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EXECUTIVE  SUMMARY 


OBJECTIVE 

This  paper  presents  a  novel  algorithm  for  estimating  the  motion  in  a  series  of  range  images. 
First,  each  range  image  is  approximated  by  applying  a  high-order  polynomial  expansion  to  local 
neighborhoods  within  the  range  image.  Then,  these  approximations  are  used  to  derive  the  transla¬ 
tion  or  displacement  estimation  from  frame  to  frame  within  the  series  of  range  images  (also  known 
as  range  image  flow).  An  iterative  method  for  computing  the  translation  is  presented. 

RESULTS 

We  evaluate  the  algorithm  on  several  synthetic  and  real-world  range  image  sequences  with 
promising  results. 

RECOMMENDATIONS 

Results  in  this  paper  are  generated  from  a  single  iteration  of  the  algorithm  in  space  and  time. 
Therefore,  our  next  step  is  to  improve  the  implementation  of  the  algorithm  so  that  using  multiple 
spatial  scales  and  past  information  improve  the  final  flow  estimation,  as  we  would  expect.  Also,  we 
will  port  the  MATLAB®  implementation  to  C++  to  improve  speed  and  meet  the  final  integration 
goal  of  the  algorithm. 
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1.  INTRODUCTION 


Estimating  motion  in  a  video  or  series  of  images  is  an  extremely  important  and  difficult  task 
in  computer  vision  and  has  numerous  applications,  such  as  autonomous  vehicle  navigation  [11]. 
Not  much  attention  has  been  given  to  the  estimation  of  motion  between  range  images,  however, 
estimating  motion  in  video  or  image  sequences,  most  commonly  referred  to  as  optical  flow,  has  a  long 
history  of  research.  The  two  most  common  and  prominent  approaches  to  optical  flow  are  known  as 
local  (or  sparse)  and  global  (or  dense)  methods.  Local  methods,  such  as  the  Lucas-Kanade  method 
[14],  estimate  the  motion  of  regions  of  interest  between  images  using  image  registration  and  warping 
techniques.  In  contrast,  global  methods,  such  as  Horn  and  Schunck’s  method  [13],  compute  a  dense 
motion  field  by  estimating  the  motion  of  each  pixel  between  images.  In  this  work,  we  are  concerned 
mostly  with  the  latter  method,  particularly  the  work  by  Farneback  [8],  which  approximates  the 
image  using  a  polynomial  expansion  of  local  patches  and  then  uses  the  polynomial  expansion  to 
estimate  the  global  displacements  between  images  for  each  pixel. 

Estimating  optical  flow  for  range  images,  also  known  as  range  flow,  is  the  main  topic  of  this 
paper.  Very  little  current  research  exists  on  this  topic,  one  of  the  few  examples  being  that  of 
Spies,  Jahne  and  Barron  [18],  however,  it  is  an  extremely  important  problem  in  a  growing  field. 
The  problem  of  range  flow  is  different  from  optical  flow  for  electro-optical  images  in  the  sense  that 
every  pixel  value  is  a  measure  of  distance  instead  of  color  or  brightness.  This  difference  makes  it 
very  difficult  to  apply  existing  and  traditional  two-dimensional  (2D)  optical  flow  methods  to  range 
flow.  For  instance,  the  brightness  constancy  constraint  used  by  many  optical  flow  methods  is  not 
valid  for  range  images.  Therefore,  our  approach  in  this  work  is  to  extend  a  well-known  global 
optical  flow  method  of  motion  estimation  based  on  polynomial  expansion  [8]  to  range  images.  We 
extend  the  method  by  using  a  high-order  polynomial  expansion  to  include  terms  in  the  z  direction 
(range  distance  to  the  sensor).  We  then  formulate  an  iterative  method  to  solve  for  displacement 
in  the  x,  y,  and  z  directions  between  range  imagery.  In  addition,  we  perform  the  calculation  at 
multiple  scales  for  robustness  and  include  displacement  estimates  from  previous  frames  to  improve 
the  overall  motion  estimation.  Promising  results  are  presented  on  several  sets  of  synthetic  and 
real-world  range  images. 
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2.  RELATED  WORK 


As  previously  mentioned,  there  is  a  vast  amount  of  research  in  optical  flow  between  color  im¬ 
ages,  usually  categorized  into  local  and  global  methods,  such  as  Lucas  and  Kanade  [14]  and  Horn 
and  Shunck  [13],  respectively.  Many  assumptions  and  constraints  have  been  introduced  in  both 
approaches  to  deal  with  noise  and  smoothness  of  the  solutions,  such  as  the  brightness  constancy 
assumption,  gradient  constancy  assumption  and  spatio-temporal  smoothness  constraints.  This  has 
led  to  a  breeding  ground  of  methods,  such  as  the  work  of  Bruhn,  Weickert,  and  Schorr  [3],  which 
attempt  to  combine  the  local  and  global  methods  to  address  the  drawbacks  and  assumptions  of 
each  individual  method.  The  most  popular  and  successful  methods  are  covered  in  more  detail  in 
two  benchmarking  papers  on  the  subject  by  McCane,  Novins,  Crannitch,  and  Galvin  [15]  and  Baker 
et  al.  [1].  The  two  papers  describe  common  databases,  procedures,  and  results  on  comparing  more 
than  20  optical  flow  methods,  with  the  Baker  et  al.  paper  being  the  most  recent  and  complete.  Of 
particular  interest  to  our  research  is  the  method  introduced  by  Farneback  in  several  papers  that  in¬ 
troduces  the  estimation  of  motion  using  the  polynomial  expansion  [5-10,  16,  17,  19].  In  Farneback’s 
work,  the  local  neighborhoods  of  each  image  are  approximated  by  a  polynomial  expansion  and  an 
analytical  solution  for  the  displacement  between  images  is  derived.  From  this  derivation,  a  robust 
algorithm  is  designed  to  compute  the  displacement,  and  thus,  motion  field,  between  two  or  more 
images  in  a  sequence.  The  method  has  proven  to  be  very  accurate  and  robust  for  2D  images  and 
has  been  included  as  a  default  algorithm  in  the  OpenCV  library  [2] . 

The  research  of  estimating  the  motion  between  range  images,  or  range  flow,  is  much  more  sparse. 
The  term  “range  image  flow”  first  appears  in  the  work  of  Gonzalez  [12],  where  he  formulates  a 
physics-based  approach  to  estimate  the  motion  of  the  range  sensor  relative  to  its  environment.  Our 
method  uses  the  same  basic  physical  model  of  the  range  sensor  as  that  used  by  Gonzalez.  One  of 
the  most  popular  and  earliest  papers  on  this  topic,  by  Spies,  Jahne,  and  Barron  [18],  notes  the 
unique  challenges  of  this  problem  and  proposes  a  basic  motion  constraint  equation  on  deformable 
surfaces.  The  constraint  solutions  are  obtained  in  a  total  least  squares  framework  and  compute  a 
dense  range  flow  field  from  sparse  solutions.  While  the  results  are  promising  for  the  range  images 
presented  in  their  paper,  the  method  is  not  directly  transferable  to  other  domains,  such  as  dense 
range  flow  in  a  moving  scene  due  to  the  large  displacements  present.  Therefore,  the  focus  of  our 
work  is  to  extend  the  polynomial  expansion  method  to  range  imagery  to  compute  dense  range  flow 
on  sequences  of  real-world  range  images. 
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3.  RANGE  IMAGE  POLYNOMIAL  EXPANSION 


In  our  formulation,  range  image  flow  uses  a  polynomial  expansion  based  approximation  of  the 
range  image.  This  approximation  is  done  using  a  set  of  quadratic  basis  functions,  applied  to  the 
range  data.  The  basis  equation  set  is  {1,  x,  y,  x2,  y2,  xy},  which  describes  the  variation  of  z,  range 
from  the  sensor,  as  you  vary  x  and  y ,  azimuth  and  elevation  with  respect  to  the  sensor.  In  addition 
to  the  basis  functions,  we  incorporate  a  notion  of  the  accuracy  or  importance  of  this  data  though 
a  certainty  matrix,  as  well  as  a  proximity-based  weight  over  the  neighborhood  in  the  form  of  an 
applicability  matrix,  as  is  done  in  Farneback’s  Ph.D.  dissertation  [8].  For  the  certainty  matrix,  a 
value  of  1  was  given  to  all  pixels  populated  with  valid  data  and  0  for  any  pixel  with  error  or  no 
data.  A  Gaussian  kernel  was  used  as  the  applicability  measure.  The  weights  of  the  basis  were 
calculated  using  Equation  (1)  as  described  by  Farneback, 

r  =  (BTWaWcB)-1BTWaWcf  (1) 

where  B  represents  the  basis  functions,  Wa  and  Wc  are  applicability  and  certainty,  taken  column¬ 
wise  and  diagonalized,  and  f  is  the  range  image  data,  taken  column- wise.  The  values  of  these  weights 
for  a  Velodyne®  and  Odetic  lidar  range  image  [4]  can  be  seen  in  Figures  1  and  2,  respectively. 
The  approximated  signal  can  then  be  reconstructed  as  Br,  reshaped  column-wise  into  its  proper 
neighborhood.  This  defines  the  value  of  each  pixel  as  a  function  of  its  local  coordinates,  x,  described 
in  Equation  (2)  as  a  quadratic  matrix  multiplication  and,  in  Equation  (3)  in  index  notation, 


/(x)  =  xrAx  +  bTx  +  c 

(2) 

/(x)  =  aijXiXj  +  biXi  +  c, 

(3) 

where  A  = 


rxy 

2 


V 


b  = 


and  c  =  r\.  While  the  matrix  form  is  more  compact  and 


easier  to  manipulate,  the  index  notation  becomes  necessary  as  the  model  moves  into  higher-order 
approximations  using  larger  degree  tensors. 
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(a)  Original  Range  Image 


(b)  1  Coefficient  Image 


(c)  y  Coefficient  Image 


(d)  x  Coefficient  Image 


(e)  y 2  Coefficient  Image 


A\ 

rf  v,  1  ■  1  '  ' 

v  .  i/  W*„< 

Wt*  *1  ' 

(f)  x 2  Coefficient  Image 


(g)  xy  Coefficient  Image 
Figure  1.  Velodyne®  polynomial  expansion. 


6 


(a)  Reflective  Image 
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(c)  1  Coefficient  Image 


(f)  y 2  Coefficient  Image 


age  (b)  Original  Range  Image 
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(d)  y  Coefficient  Image 
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(g)  x 2  Coefficient  Image 

(h)  xy  Coefficient  Image 

Figure  2.  Odetic  polynomial  expansion. 
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4.  DISPLACEMENT  CALCULATION 


The  polynomial  expansion  in  2D  images  allows  displacement  to  be  calculated  analytically  by 
looking  at  the  effects  of  a  displacement  on  the  polynomial  expansion  coefficients.  The  effects  of 
this  displacement  on  the  quadratic  polynomial  are  derived  in  Equation  (4). 

/i(x)  =  xtAix  +  biTx  +  ci 
/2(x)  =  /i(x  — d) 

=  (x  -  d)TAi(x  -  d)  +  biT(x  -  d)  +  c\  (4) 

=  xtAix  +  (bi  —  2Aid)Tx  +  dTAid  +  biTd  +  c\ 

=  xtA2x  +  b2Tx  +  c2, 

leaving  you  with  a  new  quadratic  polynomial  with  different  coefficients, 


A2  =  Ai  (5) 

b2  =  bi  —  2Aid  (6) 

c2  =  dTAid  +  biTd  +  ci.  (7) 

With  these  new  coefficients,  d  can  be  computed  from  Equation  (6): 

d  =  ~A1~1(b2-b1).  (8) 


This  method  of  displacement  calculation  was  developed  and  tested  by  Farneback  and  is  used 
as  the  starting  point  of  the  three-dimensional  (3D)  displacement  calculation. 

A  third  dimension  cannot  simply  be  added  to  d  because  with  the  current  model,  it  has  no 
meaning  as  an  input  term.  The  function  space  is  not  defined  for  any  value  of  z  under  this  model 
so  the  model  must  be  modified  to  explain  behavior  while  displacing  this  third  dimension.  One 
possible  solution  could  be  that  z  displacement  subtracts  from  the  function  value, 

"0" 

f(x-  0  )  =  /(x)  -  dz.  (9) 

_dz_ 

This  model  requires  A  and  b  be  changed  to 

"  rx2  r-f-  7i" 

A=  r~f~  V  72  (10) 

_-7i  “72  0  _ 
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(11) 


b  = 


This  model  depends  on  the  values  of  71  and  72,  which  if  viewed  as  the  coefficients  of  the  xz 
and  yz  terms,  should  be  0,  making  A  always  singular.  Otherwise  A-1  is  heavily  affected  by  a  term 
with  no  meaning  to  the  polynomial  equation. 

To  properly  capture  this  higher-dimensional  behavior,  a  higher-order  equation  is  used  to  ap¬ 
proximate  this  space.  A  linear  spreading  of  the  x  and  y  data,  Equation  (12),  as  well  as  a  constant 
increase,  Equation  (13),  was  used  to  approximate  the  behavior  of  the  data  as  you  move  along  the 
z  dimension.  This  leads  to  a  quartic  polynomial,  Equation  (14),  with  the  higher-order  terms  being 
very  sparse  tensors: 


x  =  (axz  +  l)x 

y'  =  ( OLyZ  +  1  )y 


/(x)  =  fquadi*  2)  + 

/(x)  =  axrx2X  z  +  axayr Xyxyz  +  ayry2y  z 

+  2axrx2X2z  +  (ax  +  ay)rxyxyz  +  2  ayry2y2z 
+  rx2X 2  +  rxyxy  +  axrxxz  +  ry2y 2  +  ayryyz 
+  rxx  +  ryy  +  z  +  c. 

This  can  be  combined  into  a  tensor  form  similar  to  Equation  (3): 


(12) 

(13) 

(14) 


/(x)  —  gij^iXiXjX^xi  T  hijkXiXjXk  T  aijXiXj  T  b{X{  H-  c, 


(15) 


where  the  high-order  tensors,  G  and  H,  are  sparse  and  the  quadratic  tensors  are  dense.  For  G  and 
H,  the  nonzero  terms  are 


axrx2 

90022  =  *70202  *72002  =  •••  =  - X - 

6 

^x^y^  xy  ,  \ 

9 0122  =  *70212  *72012  =  •••  =  - ^ -  (lo) 

^y^y2 

9 1122  =  *71212  =  *72112  =  •••  — - — 


^002  =  ^020  =  •••  = 


2q^7V2 


hoi2  =  h\20  =  •••  = 


(ydx  T  (Xy)T’xy 


hi  12  =  hui  =  ... 


Xy'y 2 


(17) 
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while  A  and  b  are 


rx2 

rxy 

2 

rxy 

2 

V 

CHx'T’x 

°iyry 

2 

2 

°txr  x 


2 

OLyT  y 
2 


o 


b  = 


We  now  explore  the  effects  of  displacement,  as  Farneback  did,  in  Equation  (4): 


(18) 


(19) 


f ( x)  (Jlj k I X j X j X k X I  hijfcXiXjXfc  +  (IjjXjXj  +  bj Xj  +  c 

/(x)  =  /(x-  d) 

=  9ijkl{xi  -  di)(xj  -  dj)(xk  -  dk){xi  -  di ) 

+  hijk{xi  -  di)(xj  -  dj)(xk  -  dk ) 

"F  &ij {x%  di)(Xj  dj) 

+  bi{xi  —  di)  +  c 

=  9ijklXiXjX]^Xi 

{^Qijkl^l  hijk)xixjXk 
“I-  (6 Qijkl  Shijkdk  +  aij)xiXj 

{^dijkldjdkdi  Shijkdjdk  e2aijdj  bi)xi 
(dijkldidj dfcdi  hijkdidjdk  H-  aijdidj  b{di  -\-  c 

—  9ijklX%X jXj^Xi  ~\~  hijkXiXjXk  +  CLijXiXj  b{X{  -\-  c, 

leaving  us  with  a  new  quartic  polynomial  with  the  following  coefficients: 


ii 

to. 

(21) 

1 

-cT 

II 

(22) 

dij  CLij  3hijkdj~  +  6gijkidkdi 

(23) 

hi  =  hi  ^ttijdj  H-  ^hijfcdjdk  Agij^idjd^di 

(24) 

c  —  c  bidi  ctijdidj  hijkdidjdk  H-  Qijkldidjd^di.  (25) 

Now  we  have  a  linear  equation  with  H  that  could  be  used  as  an  analytical  solution  to  d  but 
with  our  higher-order  tensors  being  mostly  sparse  and  mainly  composed  of  the  model  coefficients 
as  opposed  to  the  expansion  coefficients,  we  look  to  our  lower  order  dense  tensors  to  solve  for  d 
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through  numerical  optimization.  Using  the  symmetries  and  sparsities  of  the  matrices,  we  are  left 
with  the  following  coefficients  in  our  dense  tensors.  In  this  case,  62  =  1  because  it  is  created  purely 
by  the  model  of  z  motion  and  has  information  about  the  next  image’s  polynomial  expansion: 


bo  =  bo  —  2(aoodx  +  ao\dv  +  ao2dz) 

+  6(h002dxdz  +  hoi2dydz) 

—  12(50022  dxd2z  +  goi22dyd2z) 
bi=bi  —  2(a0idx  +  audy  +  audz) 

+  §{hoi2dxdz  +  hn2dydz)) 

—  12(goi22dxdl  +  gu22dyd2z) 
b2  =  b2-  2 (a02c4  +  audy) 

+  3(hoo2  d^x  +  h\ud2  +  2houdxdy) 

—  12(50022^^2  +  9ii22d2dz  +  2g0u2dxdydz) 

=  1 


c  =  c  —  bodx  +  bidy  +  b2dz 

+  (uoo^  ~b  and2  +  2ao\dxdy  +  2cio2dxdz  +  2a,ndydz) 
-  2>(h0o2d2xdz  +  hiud2ydz  +  2h012dxdydz) 

+  6(50022^  +  gii22dyd2  +  2goi22dxdyd2z) . 


(27) 


To  calculate  d,  we  optimize  the  difference  between  the  observed  polynomial  coefficients  of  the 
next  image  and  the  coefficients  derived  through  displacing  the  coefficients  of  the  first  image  using 
Equations  (26)  and  (27).  This  optimization  is  done  using  the  non-linear  least  squares  version  of 
Newton’s  method,  the  Gauss — Newton  algorithm,  on  Equations  (28)  and  (29).  This  optimization 
technique  requires  the  Jacobian  matrices  of  each  coefficient  with  respect  to  d,  calculated  using 
the  partial  derivative  in  Equations  (30)  and  (31).  Again,  b2  contains  no  polynomial  expansion 
information,  so  its  components  of  the  Jacobian  are  0. 

mm£(&<1|(d)-i,f))2  (28) 


min(c(’1^(d)  —  c^)2 

d 


(29) 
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db0 

ddx 

dbi 

ddx 

dh 

ddx 

db0 

ddy 

dbi 

Wy 

db2 


ddy 

dbo 

ddz 

d'h 

ddz 

db2 

ddz 


— 2aoo  +  6W*  —  3^0022^ 

— 2aoi  +  6hoi2dz  —  ‘igovndi 

— 2ao2  +  6(/ioo2^x  +  hoi2dy)  —  6(goo22dxdz  +  goi22dydz)  —  0 
— 2aoi  +  6/1012  dz  —  3goi22d2z 
— 2an  +  6/iii2^4  —  3^1122^ 

— 2ai2  +  6(4oi24  +  hn2dy)  —  Q(goi22dxdz  +  gn22dydz)  —  0 
— 2ao2  +  6(W*  +  howdy)  —  Q(go(mdxdz  +  goi22dydz) 

— 2ai2  +  6(4oi24  +  h\i2dy)  —  6(goi22dxdz  +  gn22dydz) 

—  3(2^0122^^2/  +  g0022dl  +  31122^)  =  0 


(30) 


dc 


dc 

ddy 

dc 

ddz 


—  bo  +  2(aoo<4:  +  “1“  a02 dz) 

-  6(4oo2 dxdz  +  hoi2dydz )  +  l2(go(y22dxd2  +  goi22dyd2z) 

— 6i  +  2(uoi<4;  3“  ®ii^3/  T  a  +  12<4) 

-  6(40i2^4  +  hn2dydz )  +  12(goi22dxdz  +  gn22dyd2z) 

—  1  +  2(ao2dx  +  a\2dy)  —  3(hoo2d2  +  h\i2d 2  +  2ho\2dxdy ) 
+  12(50022^^  +  9ii22d2dz  +  2g0i22dxdydz ) 


(31) 


This  leaves  two  remaining  routes  to  solving  for  d,  through  the  changes  in  b  or  the  changes  in  c, 
shown  in  Equations  (32)  and  (33),  respectively,  each  with  its  own  advantage.  The  solution  based 
on  b  tends  to  produce  more  accurate  results  for  dx  and  dy  because  it  captures  the  motion  of  the 
quadratic  components  seen  mostly  in  edges.  This  solution  does  not  tend  to  be  as  accurate  for  the 
dz  component  because  dz  is  detected  through  the  spreading  of  the  quadratic  x  and  y  components. 
This  effect  can  be  small  at  long  distances,  where  as  dz  has  a  large  effect  on  c.  However,  the  effects 
on  c  do  not  track  the  quadratic  components  of  dx  and  dy  as  directly  as  b.  To  maximize  the  accuracy 
of  d,  a  weighted  average  of  db  and  dc,  shown  in  Equation  (34),  with  [ixy  relatively  high  and  ftz 
low,  where  as  7 xy  is  relatively  low  and  7,  high. 


Adb  =  (JbTJb)"1JbT(b(1)(d)  -  b(2))  (32) 

Adc  =  (JcTJc)-1JcT(c«(d)  -  c(2>)  (33) 
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fixy 

Ixy 

Ad  = 

fixy 

Adb  + 

Ixy 

Jz_ 

_lz  _ 

(34) 


Similar  to  Farneback,  this  algorithm  estimates  the  displacement  over  a  neighborhood  /  around 
x  as  opposed  to  a  single  point,  minimizing  the  following  equations: 


^  u.(Ax)||(Si1)(x  +  Ax  -  d)  -  6f(x  +  Ax))  ||2 

Axel 

w(Ax)||(c(1)(x  +  Ax  -  d)  -  C(2>(x  +  Ax))||2, 

Axel 

where  w{ Ax)  is  the  neighborhood  weighting  function  and  the  minimum  steps  are 


Adb 


u>(Ax)JbTJb 


Y  w(Ax)JbT(b(1)(d)  -b(2)) 

Axel 


Adc 


ry(Ax)JcT(c^1^(d)  —  cV*). 

Axel 


(35) 

(36) 


(37) 


(38) 
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5.  RANGE  IMAGE  FLOW  EXPERIMENTS 

This  algorithm  has  been  tested  using  data  from  a  Velodyne®  HDL-64E.  This  sensor  is  a  360° 
field  of  view  3D  lidar  with  64  vertically  mounted  lasers  on  a  spinning  head.  The  lasers  have  a 
maximum  range  of  50  m  and  an  accuracy  of  2  cm.  It  is  capable  of  spinning  at  5  to  15  Hz,  generating 
over  1.333  million  points  per  second.  For  our  tests,  the  sensor  was  set  to  10  Hz,  generating  a 
horizontal  resolution  of  1800  returns  per  rotation.  The  lidar  returns  are  assembled  into  a  64  x  1800 
range  image.  The  tests  were  done  using  a  Ford  F-150  with  the  Velodyne®  mounted  to  the  roof.  The 
vehicle  and  mounting  hardware  are  visible  in  the  lidar  scans,  so  all  pixels  within  a  threshold  have 
been  marked  with  a  certainty  value  of  0,  causing  these  values  to  have  no  effect  on  the  polynomial 
expansion  or  the  flow  calculations. 

The  data  set  shown  in  Figures  5,  6,  7,  8,  and  9  was  collected  along  Woodward  Road  at  SSC 
Pacific  on  a  foggy  day.  The  images  contain  mostly  data  from  the  bushes  and  other  vegetation 
surrounding  the  road,  but  do  eventually  show  a  parking  lot,  seen  in  Figure  9.  No  ground  truth  is 
currently  available  for  the  flow  fields  so  only  qualitative  analysis  could  be  done. 

Most  of  the  flow  field  appears  accurate,  though  some  regions  still  contain  peculiar  behavior. 
Certain  regions  of  the  flow  image  act  as  sources  or  sinks  to  the  flow  fields,  such  as  in  Figure  6, 
where  in  the  middle  right  portion  of  the  image  the  flow  field  moves  away  from  a  central  point.  The 
boundary  of  the  scan  also  seems  to  be  less  accurate  than  the  central  areas.  Areas  such  as  that 
shown  in  Figure  8  are  relatively  uniform,  making  flow  calculations  more  difficult.  Despite  this, 
the  algorithm  performed  reasonably  well,  though  it  does  have  regions  with  sporadic  flow  behavior. 
Further  work  and  testing  will  be  done  to  this  algorithm  to  improve  its  performance  and  a  ground 
truth  test  will  be  run  to  quantitatively  evaluate  the  performance  of  this  algorithm. 


14 


6.  POLYNOMIAL  CLASSIFICATION  AND 
PLANAR  REGION  DETECTION 


Using  Odetic  Laser  Range  Finder  range  images  [4]  of  a  series  of  geometric  objects,  we  were 
able  to  test  the  polynomial  coefficients  associated  with  3D  features.  The  coefficients  associated 
with  specific  3D  features  are  shown  in  Table  1.  The  primary  features  of  interest  are  planes,  edges, 
and  corners.  Figure  3  shows  the  comparison  between  the  real  image  patches  (3a,  3c,  and  3e) 
associated  with  these  features  with  their  respective  polynomial  expansions  (3b,  3d,  and  3f,  respec¬ 
tively).  Planar  features,  such  as  the  one  shown  in  Figure  3a,  have  quadratic  components  close 
to  zero  and  linear  components  dependent  on  the  normal  of  the  plane.  This  planar  feature  has 
an  expansion,  {rq,  ry,  rXl  ry2,  rq2,  rxy},  of  {68.0,0.19,-0.11,0.00,0.02,-0.01}.  The  polynomial  ex¬ 
pansion  shows  the  very  low  quadratic  components  indicative  of  planar  regions.  Edge  features, 
however,  have  quadratic  features  depending  on  the  direction  of  the  edge  and  linear  features  de¬ 
pendent  on  the  location  of  the  edge.  In  Figure  3c,  there  is  a  horizontal  edge  with  an  expansion  of 
{85.7,-11.19,0.26,2.95,0.01,-0.10}.  In  this  case,  the  expansion  shows  the  relatively  high  x  and 
x 2  components  found  in  horizontal  edges.  Finally,  corner  features  have  nonzero  zero  values  for  all 
of  their  linear  and  quadratic  components.  The  expansion  of  Figure  3c  shows  an  example  of  this, 
{137.1,  —19.28,  24.1,  2.13,  2.64,  7.53},  showing  all  nonzero  linear  and  quadratic  terms. 

Table  1.  Polynomial  Expansion  of  3D  Features. 


Polynomial  Coefficients 


ry 

rx 

V 

rx 2 

rxy 

Feature  Type 

any 

any 

ft  0 

« 0 

«  0 

Plane 

ft  0 

ft  0 

96  0 

0 

ft  0 

Vert.  Linear 

96  0 

ft  0 

96  0 

« 0 

«  0 

Vert.  Edge 

ft!  0 

ft  0 

ft  0 

96  0 

ft  0 

Horz.  Linear 

ft  0 

96  0 

ft  0 

96  0 

«  0 

Horz.  Edge 

ft  0 

ft  0 

any 

0 

96  0 

Diag.  Linear 

96  0 

ft  0 

any 

0 

96  0 

Diag.  Edge 

ft  0 

ft  0 

ft  0 

any 

96  0 

Diag.  Linear 

ft  0 

96  0 

ft  0 

any 

96  0 

Diag.  Edge 

ft  0 

ft  0 

96  0 

96  0 

any 

Spherical 

96  0 

96  0 

96  0 

96  0 

any 

Corner 

This  leads  to  the  development  of  a  planar  segmentation  algorithm.  Using  the  polynomial 
expansion  of  the  range  image,  the  planar  regions  can  be  differentiated  from  the  non-planar  regions 
by  the  quadratic  terms  of  the  expansion.  From  there,  a  region-growing,  flood-fill  technique  is  used 
to  determine  which  regions  belong  to  the  same  plane.  In  this  case,  the  pixels  at  the  frontier  of  the 
region  are  checked  for  low  quadratic  components  and  if  they  are  found  to  be  planar,  then  their 
linear  components  are  compared  to  the  frontier’s  linear  components.  If  the  linear  components  are 
sufficiently  different,  the  pixel  becomes  the  seed  to  a  new  region;  otherwise,  it  is  added  to  the 
frontier  and  the  region  continues  to  grow.  This  method  was  tested  on  the  Odetic  range  images  and 
an  example  is  shown  in  Figure  4.  The  floor  is  correctly  segmented  as  a  single  plane  and  the  wall 
is  correctly  segmented  into  two  planes,  with  a  break  caused  by  a  raised  wire  covering.  The  box 
structure  is  also  correctly  segmented. 
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Figure  3.  Outlined  planar  (a),  edge  (c)  and  corner  (e)  features  with  their  respective  polynomial  expansions  (b,  d 
and  f,  respectively). 


Figure  4.  Planar  region  detection  of  Odetic  range  image.  Non-connected  colors  are  independent  planar  regions. 
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Figure  5.  Range  image  flow  on  Velodyne®  scan  100. 
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Figure  6.  Range  image  flow  on  Velodyne®  scan  200. 
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Figure  7.  Range  image  flow  on  Velodyne®  scan  204. 
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Figure  8.  Range  image  flow  on  Velodyne®  scan  400. 
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Figure  9.  Range  image  flow  on  Velodyne®  scan  1000. 
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7.  CONCLUSION 


The  Polynomial  Expansion-based  method  is  effective  in  calculating  the  flow  of  range  data, 
though  further  testing  should  be  done.  In  addition,  this  method  can  detect  simple  3D  features 
and  calculate  planar  segmentation  on  range  image  data.  In  future  work,  quantitative  testing  will 
be  done  to  validate  the  qualitative  results,  using  ground  truth  data  to  run  large  scale  quantitative 
tests.  The  feature  detection  and  planar  segmentation  can  be  incorporated  into  the  flow  calculation. 
Certain  features  have  more  stable  information  about  their  motion  and  as  such  should  take  on  greater 
weights  in  the  flow  calculations.  Planer  regions,  for  example,  contain  very  little  information  about 
the  data’s  motion,  but  the  edge  and  corner  features  contribute  a  large  amount  of  flow  information. 
Using  the  planar  segmentation,  the  flows  calculated  at  the  corners  and  edges  of  planar  regions  could 
be  interpolated  across  the  region.  In  addition  to  the  range  image  work,  future  work  will  focus  on 
calculating  the  rotational  components  of  the  flow  vectors  with  respect  to  a  single  image  location, 
useful  in  image  stabilization  and  unmanned  aerial  vehicle  orientation  determination. 
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